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Abstract 

A new Monte Carlo event generator for wide angle Bhabha scattering and muon 
pair production in e + e~ collisions is described. The program includes complete 
one-loop electroweak corrections, and QED radiative corrections. The 0(a) QED 
correction uses the exact matrix element. Higher order QED corrections are included 
in an improved soft photon approximation with exponentiation of initial state ra- 
diation. Events are generated in the full phase space of the final state including 
explicit mass effects in the region of collinear mass singularities. The program is 
intended for centre of mass energies around and above the Z peak and for Bhabha 
scattering at angles greater than 10° . 
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PROGRAM SUMMARY 

Title of the program: BHAGENE3 

Computer: IBM 3090 

Operating system: VM/CMS 

Programming language used: FORTRAN 77 

High speed storage required: 678k words 

No. of bits in a word: 32 

Peripherals used: Line printer 

Number of cards in combined program and test deck: about 6900 

Keywords: Radiative corrections, Monte Carlo simulation, Wide angle Bhabha scat- 
tering, Muon pair production, Photons, Quantum electrodynamics, Electroweak theory 

Nature of physical problem: Calculation of all 1-loop electroweak corrections, and of 
QED corrections up to 0(a 3 ) for wide angle Bhabha scattering and muon pair production 
in the vicinity of the Z peak. 

Method of solution: The 1-loop electroweak virtual corrections are calculated analyt- 
ically. The 0(a) QED radiative corrections use the exact matrix element. The higher 
order QED corrections are given by an improved soft photon approximation. Event con- 
figurations containing real photons are first generated according to an approximate model, 
and are then reweighed according to the theoretical distributions. The weight throwing 
technique is used to produce unit weight events in the full final state phase space. 

Restrictions on the complexity of the problem: Some terms in vnijE (I = e, /i) are 
neglected, so that fermion threshold effects are not properly taken into account. The 
approximation used for 0(a 2 ) and higher order real photon radiation may be unreliable 
in events with multiple hard photon radiation. For cross section accuracy at the % level 
it is not recommended to use the program for Bhabha scattering angles less than 10°. 

Typical running times: For Bhabha events: initialisation about 140 sec, generation 709 
unit weight events per second. For muon pairs the initialisation time is shorter by a factor 
of about 3 and the event generation rate is roughly doubled. The time unit is an IBM3090 
CPU second. 

Unusual features of the program: Extensive use is made of one and two dimensional 
look-up tables for fast and flexible generation of Monte-Carlo variables. The space re- 
quired for these tables means that the fast memory requirements may be greater than 
in other comparable programs. These look-up tables as well as average event weights 
are created in the initialisation phase of the program which is, in consequence, relatively 
time-consuming. 
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1 Introduction 

Data on the leptonic decays of the Z taken since 1990 at LEP and SLC have made 
possible many precise tests of the Standard Model (SM) of weak and electromagnetic 
interactions 0. These electroweak analyses require Monte Carlo event generators, in- 
corporating both the SM predictions and the numerically important pure QED radiative 
corrections, in order to properly take into account detector acceptance, efficiency and 
resolution effects, for arbitrary experimental cuts. 

For fermion pair production (excluding e + e~), the most widely used generator is KO- 
RALZ 0. This incorporates one- loop electroweak corrections using the SM, multiple 
bremsstrahlung photons in the initial state and a single bremsstrahlung photon in the 
final state as well as a complete simulation of r decays, including polarisation effects and 
QED radiative corrections associated with the r decay products. 

For wide angle Bhabha scattering two generators have been written prior to that, 
BHAGENE3, described in the present paper. The first, BABAMC @, g includes one- 
loop SM electroweak corrections and QED corrections (real and virtual) to 0(a). The 
second, BHAGEN || includes higher order QED corrections (by exponentiation) in both 
initial and final states, but is limited to 'quasi-elastic' configurations where the photons 
are soft, or almost collinear with the radiating lepton. This limitation is a severe dis- 
advantage when comparing with actual experimental data, since in general very loose 
cuts (typically none at all on the lepton-photon angles) are made in order to minimise 
the size of the QED radiative corrections. More recently a new Monte Carlo generator 
UNIBAB || has been written, which is a general purpose program in many ways com- 
parable to BHAGENE3. Electroweak corrections are included and arbitrary numbers of 
initial and final state photons generated using a LLA structure function formalism. Un- 
like BHAGENE3 however, 0(a) initial/final interference effects are not included in the 
current version. 

Bhabha generators also exist |7], |] for the small angle region of lepton scattering 
angle Q[ (dominated by the t-channel photon exchange diagram) typical of the angular 
acceptance of the luminosity monitors of the LEP/SLC experiments. 
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The generator described in this paper, BHAGENE3, is adapted to the description of 
wide angle Bhabha scattering (10° < Q\ < 170°), and of /x-pair production, in the region of 
the Z peak, r-decays are not incorporated. Electroweak corrections in the SM, including 
the most important two- loop effects, are implemented as described in Ref.|J. Multiple 
hard photon generation, ( n 7 < 2(3) for initial (final) states) is also included. No restric- 
tions are necessary on the kinematical configuration of the generated events. However, 
the approximate nature of the hard photon generation algorithm should be borne in mind 



when multiple, very hard photon configurations are considered [l0|. A brief description 



of the program has been given previously, and comparisons of cross-sections and charge 



asymmetries made [ 10 1 with the analytical or semi- analytical programs: ZFITTER [] I 



ALIBABA [0 and TOPAZ0 @. 

The approach adopted for the QED radiative corrections is to treat the 0(a) correc- 
tion exactly (in particular the relevant fermion mass terms are systematically included 
throughout, so that collinear photon radiation is correctly generated) as in BABAMC 0, |4]] 
and MUSTRAAL ||14 |. Higher order real photon radiation is treated, not by using exact 



0(a 2 ) amplitudes [TS|, [T(| O, but by an improved soft photon approximation. This may 
be justified |IJ by the relatively small size of the 0(a 2 ),0(a 3 ), ... hard photon corrections 
as compared to that at 0(a). As in Ref. [3,4,14] the final state phase space is divided into 
two regions. A cut (Eq ) is applied on the energy of each photon. For E 1 < Eq the photon 
energy and angular distributions are integrated over, so as to yield a Virtual, Soft (VS) 
corrected cross-section. The corresponding generated events have a 'Born topology'; the 
outgoing leptons being exactly back-to-back with energy: Ei = E beam = E. The gener- 
ation of such V,S events is described in Section 3 below. For E 1 > Eq a dilepton event 
with 1,2,3 hard photons is generated over the full available phase-space. Denoting the 
number of initial/final state photons by n^/n^, the different possibilities are: 

n\/n* = 1/0, 0/1, 2/0, 1/1, 0/2, 0/3 
This 'hard photon' generation is described in Section 4 below. 



A complete description of the program, its structure and parameters, may be found 
below in Section 5. In this Introduction only a brief account is given of the most important 
techniques employed (see also Ref. [10]). The first step in the execution of the program 
is an initialisation phase where all electroweak parameters are calculated within the SM 
from standard input parameters such as the fine structure constant, a, the Fermi constant 
G, and the masses of the Z (M z ), the top quark (M t ) and the Higgs boson (M H ). Pa- 
rameters needed for the QED radiative corrections, taking into account the lepton flavour 
(e,/z) selected and the beam energy E are also calculated at this stage. Next the proba- 
bilities for VS, 'Initial State hard' (IS) and 'Final State hard' (FS) events are calculated. 
For the VS events this is done by angular integration of the differential cross-section, 
including weak and electromagnetic radiative corrections. For IS, FS events approximate 
factorisable cross-sections are used that may be readily integrated analytically and/or 
numerically, to give the corresponding probabilities. Also, during the initialisation phase, 
a number of one or two dimensional Look Up Tables (LUT) for the total photon energy, 
photon energy splitting fractions, and lepton scattering angle, are produced. For each 
variable the Integrated Probability Distribution (IPD) is calculated by analytical and/or 
numerical integration. The IPD is then inverted by linear interpolation (see Appendix 
C) to yield the LUT. The average weights Wj, W p of the IS, IF events are also found, 
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during initialisation, by re-weighting events according to the 'exact' hard cross-section, 
to be described below in Section 4. 

The event generation phase of the program execution is now entered. A VS , IS or 
FS event is first chosen. For IS, FS events or are then chosen according to a 
Poisson distribution, and the appropriate sub-generator is entered. The extensive use 
of LUT rather than a Weight Rejection Procedure (WRP) as in Ref.[4,13] gives a very 
fast and efficient event generation algorithm. Arbitrary user defined cuts may be applied 
during the event generation phase. Four-vectors are written out for each unit weight event 
which is chosen by using a WRP. In the last stage of execution the 'exact' cross section 
corresponding to the user-supplied cuts is is calculated and printed out as well as the 
user-defined histograms or other distributions that are up-dated during the generation 
phase. 

Finally in this Introduction a few remarks on the limitations of the program and 
some recommended restrictions on its use. Only cross-sections summed over final and 
averaged over initial states of the lepton and photon polarisations are calculated. Only 
leptonic (not quark anti-quark) final states may be generated. As in Ref. [3,4,14] not all 
terms ~ mi/E are included, so that the near-threshold region of lepton pair production 
is not properly described. As an approximate model is used to generate hard photons 
at 0(a 2 ) and higher, based on an improved soft photon approximation, cross-sections 
for configurations with multiple hard photons should be regarded with caution. Some 
comparisons with exact 0(a 2 ) calculations are given in Ref. [10]. The precision of the 
LUT's used for the lepton scattering angle requires that 10° < Q\ < 170° in Bhabha 
scattering, if a cross-section accuracy of ~ 1% is required. 



2 Weak Virtual Corrections 



In the initialisation, the W mass is determined iteratively from a, G^, M z in subroutine 
SETCON: 



Mi 



w 



M, 



\ v^Mf (1 - Ar) 
where Ar is calculated in subroutine SEARC1. 



1 - 



1 - 



A-KQL 



(2.1) 



The improved Born approximation for the Bhabha differential cross section consists 
of the sum of contributions with t channel and s channel exchange of the photon and Z 
boson and their interferences: 



da wc na 2 



dc 



2s 



J2J2 T «( A ) A = 1 , 1 Z,Z, a = s,st,t 



(2.2) 



A a 



Here s and t are the usual Mandelstam variables and c = cos Q\ where Q\ is the lepton 
scattering angle. In Bhabha scattering the final state helicities are not measurable, while 
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the beams may be longitudinally polarized. We introduce the notation: 



Ai 
A 2 
A 3 



1 - A+A_ 
A+- A_ 
1 + A+A_ 



(2.3) 
(2.4) 
(2.5) 



where A+(_) is the degree of longitudinal polarisation of the positron (electron) beam. 
The s channel cross section contributions are: 



T s (l) 
T S {Z) 



Ai (l + c- 



\Fa(s)\ z 

2^e{p(s) X (s)F* A (s) [[Aitfee(s) + A 2 i;e(s)](l + c 2 ) 
+ [(Ai + A 2 K(s)]2c]} 

|p(^)| 2 |[Ai(l + 2K( S )| 2 +| Wee ( s )| 2 ) 

+ 2A 2 3?e (v e (s)[l + v ee (sT])](l + c 2 ) 

+ 23?e [Ai (\v e (s)\ 2 + v ee (sj) + \ 2 v e (s) (1 + v* ee (s)) 



(2.6) 
(2.7) 



2c 



Here, we use the following abbreviations: 

x(«) = 



s - M 2 z + isT z /M z 
M 2 Z 



V2 8 



ira 



- 4s w K e (s) 
-1 + 2f e (s) + 16s^K ee ( 



(2.8) 



(2.9) 

(2.10) 

(2.11) 
(2.12) 



In our expression ( |2.2|) for da wc / dc the axial couplings are equal to unity. The effective 
couplings v e (s),v ee (s) and the weak form factor p(s) are complex valued. They contain 
virtual weak corrections which are called from the electroweak library BHASHA which is 
part of the program package described here. It has been derived from the electroweak 



library DIZET |18[ which is used in the package ZFITTER |TT| and is intended for the 
description of s channel fermion pair production in the region of the Z resonance. The 
complete theoretical description of the renormalisation procedure adopted in the unitary 
gauge and related topics may be found in [19|, the formulae for the Z width in jMj and 
s channel scattering in [IT], and those for Bhabha scattering in M. The complex valued 
s channel corrections from the fermionic vacuum polarisation are contained in Fa(s), 



Fa(s) 



Q( 



\s\) 



a 



(2.13) 



and are explained below. 



We now describe the t channel contributions: 

1 



T t ( 7 ) = F A (ty 



2A1 (I^ +8A3 (I^ 



(2.14) 
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T t ( 7 Z) = 2p(t) X (t)F A (t) 



x { 2 [A x (1 + v ee (t)) - \ 2 v e (t)\ £ + % - 8A 3 (1 - v ee (t)) 



(1-c) 2 ° v ^"(l-c) 



2 



(2.15) 



T t (Z) = [p(t) X (t)] 2 ^2 

+A\ 2 v e (t) (1 + v ee (t)) 



Ax (l+4t; e (t) 2 + 2t; ee (t)+t; ee (t) 2 ) 
1 + c) 2 



8\ 3 (l-v ee (t))- - (2.16) 



;i- c ) 2 ov eev yy (1-c) 5 

The following additional abbreviations are used: 

X (t) = k l — v (2.17) 

t = ~(l-c) (2.18) 

The form factors F^(t) , v e (t) , v ee (t) , p(t) are real valued. For wide angle Bhabha scattering, 
the values of \t\ may become substantially smaller than s. Nevertheless, the form factors 
do not vary much since they depend only logarithmically on the scale. There is one 
exception to this. Bremsstrahlung diagrams with initial state radiation and t channel 
photon exchange may yield substantial cross section contributions. There, the effective 
value \t'\ of \t\ may become extremely small and the different value of the running QED 
coupling has to be taken into account properly: 

FAt) = (2.19) 
a 

For this reason the running alpha correction is included not only in da wc / dc but also in 
the cross sections for final states with hard photons (see below). 

Finally, we describe in this section the contributions from the 7Z interference: 

T st ( 7 ) = -2^e[F* A (s)F A (t)]X 1 { -^^ (2.20) 

(1-c) 

T^Z) = -2^e{x(t)p(t)FX(s)[\ 1 (l + vUt))~2\ 2 vM + (t^s)}^^ 

(2.21) 

T st (Z) = -2^e{ X (s)p(s)x(t)p(t)[X 1 ([l + v ee (s)}[l + v ee (t)} 

+4v e (s)v e (t)) - A 2 (v e (s)[l + v ee (t)\ + v e (t)[l + v ee (s)])]} (1 + C) 

(1 - c) 

(2.22) 

The running of the QED coupling a(Q 2 ) is taken into account as follows: 

««») = ^ (2.23) 
Aa = Aai + Aa udcsb + Aa t (2.24) 

where we use the convention that in the s channel it is Q 2 = — s, and in the t channel 
Q 2 = \t\. The Aa is calculated in function XFOTF1. 
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For leptons: 



AF f {Q 2 ) 



E Q 2 f N f AF f (Q 2 

f=e,fi,T 



a: 

71 



+ 



4 m 



9 3Q 2 3 



2m^ 



In 



i7r6{-Q 2 - 4m^ 



\ 



1 + 



4m 2 



(2.25) 

>]) 

(2.26) 
(2.27) 



and their colour factor Ay and charge Qf are unity. In the weak library, the AFf is 
calculated by function XI3, AFf = 2 XI3. 

The contribution from the light quarks has been parametrised in two different ways. 
Either it is calculated with ( |2.25| ) (with flag setting NPAR(2)=2); in this case, we use 
effective quark masses [ffiH : m u = m d = 0.041, m s = 0.15, m c = 1.5, m b = 4.5 GeV. 
The other, preferred approach (with NPAR(2)=3) uses a parametrisation of the hadronic 
vacuum polarisation which is contained in function XADRQQ. 



The t quark corrections are: 

Aa t = Q 2 t N t AF t {Q 2 ) + Aa 2loop > aa * 



(2.28) 



where the latter term contains higher order corrections and is calculated from functions 
ALQCDS, ALQCD: 



Aa 21oop '' 



aa, 2 mr \ VF 2 45 



3tt ; 



Q 



4 



(2.29) 



IlY F (Q 2 ) is a two loop self energy function [[E], 23]. Setting the corresponding flag 
NPAR(3) to zero (default value) results in the neglect of this very small correction. For 
NPAR(3) = 1,2 approximate, exact calculations according to Eqn.(2.29) are made. The 
first one is not a really good approximation, the second is very time consuming. 

For the function AFf(Q 2 ) one may derive the following two approximations. For light 
fermions with m 2 ? << \Q 2 \, the following approximate formula is valid: 



a 



07T 



hi 



\Q 2 



°--ik9{-Q 2 ) 



(2.30) 



While, heavy fermions with mj >> \Q | practically decouple: 



2 ' 



At LEP 1, the effective QED coupling may be treated as a constant in the s chan- 



nel [24|: 



„ 137.036 



7 



The form factors K e , K ee , p are calculated in subroutine ROKAP(. . ., s, t, u,. . ., QE, 
QF, . . ., XFF,. . .). XFF is a vector of 4 functions! which depend on s,t,u and the charges 
of the initial and final state fermions (here minus one): 



and 



u 



-t 



(2.33) 





= XFF(1; 


u, 


—s 


t 


-1, 


-1) 


v e {s) 


= XFF(2; 


u, 


—s 


t 


-1, 


-1) 


Vee(s) 


= XFF(4; 


u, 


—s 


t 


-1, 


-1) 


P(t) 


= XFF(1; 




-t, 


u 


-1, 


-1) 


V e (t) 


= XFF(2; 




-t, 


u 


-1, 


-1) 


V ee (t) 


= XFF(4; 


s, 


-t, 


u 


-1, 


-1) 



(2.34) 



(2.35) 

These corrections are switched on and off with flag NPAR(l); the order o?m\ contributions 
to them with NPAR(6), the aa s m 2 corrections with flag NPAR(3), and the ZZ and WW 
box terms with NPAR(4). The latter are negligibly small at LEP la. They introduce 
in the weak corrections to the s channel a dependence on the scattering angle. In the 
t channel, correspondingly, the weak corrections will depend not only on the scattering 
angle, but also on s. 



The total Z width [2~0| is used in ( |2.9p and is calculated in subroutine ZWRATE: 



I2ix 



QCD 
f 



\ 



Amj 



3 a 



1 + -^Qf)pf 



X 



1 + 



2m 2 / 
Ml, 



l-A K z f s 2 w \Q f \+8( K z f ) 2 s 4 w Q} 



3m 2 / 
M 2 



(2.36) 



The form factors p z , k z are not identical with the process dependent form factors for 
the scattering process although the leading terms of the latter, if calculated at the Z 
peak, yield a good approximation of the first ones. The overall factor Nc(f) is equal to 
1 for leptons and 3 for quarks and R^ CD describes final state QCD corrections in case 
of quarks: Rf D = 1 + a s /n + . . .. The numerical value has to be chosen in accordance 
with the scheme and order of the QCD corrections in mind and is set by the user with 
the parameters XPAR(IO) for the b quark and XPAR(9) for the other, light quarks. More 
details on the Z decay rate may be found in HTT|, [25f and references therein. 



3 Virtual and Soft Photonic corrections 



Three distinct contributions to the virtual and soft radiative correction may be dis- 
tinguished: 

3 The third of these four functions, Vf, is for Bhabha scattering equal to the second one, v e . 
4 This statement depends to a certain extent on the calculational scheme chosen. 
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1) Purely weak virtual corrections as described in the preceeding Section. These 
include all self-energy and vertex corrections involving W, Z and H (Higgs) bosons, as 
well as ZZ and WW box diagrams. 

2) Electromagnetic virtual corrections: vertex corrections involving only photons, 77 
and 7Z box diagrams and the effect of both leptonic and hadronic Vacuum Polarisation 
Insertions (VPI) in all off-shell photon propagators. 

3) Corrections due to real soft photons. 

For 1) BHAGENE3 uses the 'Dubna-Zeuthen' (DZ) renormalisation scheme described 
in Refs. and [9]. The weak virtual corrections modify the vector and axial- vector 
coupling constants appearing in the s, t channel Z-exchange diagrams according to Eqns. 
3.15-3.18 of Ref.[9]. For 2), the 0(a) vertex and 77 and 7Z box contributions are taken 
from Ref.[3], and the VPI contribution from Ref.[27]. Vertex corrections to 0(a 2 ) and the 
corresponding exponentiated soft photon correction, 3), (the integral over the angles and 
energies of all photons with total energy below a fixed cut-off: y = E] ot / E < y ) are given 
by Ref. [28,29]. The differential cross-section, corrected for weak loop effects 1), and also 
the VPI part of 2), is denoted by d wc a /dc where c = cos 61 , and Oi is the lepton scattering 
angle. Including the soft real photon, the 0(a) virtual, and the leading log 0(a 2 ) virtual 
corrections as well as the 77, 7Z box diagrams, gives the following expression for the 
Virtual, Soft differential cross- sect ion: 



da 



vs 



dc 



dc 



dc 



exp(/5 e In — 
Vo 



1 



+ 



da 



EC 



dc 



Cl + {f3 mt + l3 f )\n 



Vo 



N J EC 
_|_ i fiEBO.X 

i=l 



dc 



(3.1) 



where 



Pe 

Pf 

Pint 

m f 



2a 

7T 

2a 

TT 



"'4 

mi 



In 



- 1 



m, 



4a , t 
— In — 

7T u 

1 + 2! 

7T 



1 + 



a 



TT 



3 



mi 



ln(- 



m 



f 



2 

3 

7T 2 

2 

3 



mass of final state leptons 



' TT 



mi 



12 Pe 



The differential cross section da EC / dc is corrected for the 'running' of a in the s and 
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t-channel photon exchange diagrams by the replacement: 



a — > a(u) = % , . (3.2) 

where a is the on-shell QED fine structure constant (a -1 = 137.036...) and Il 7 (n) = —Aa 
is the photon proper self energy function, parameterised according to Ref. [27]. The most 
important virtual weak corrections are also included in da EC / dc via the replacements 

where is the Fermi constant, determined from the muon lifetime, and sw = sin 6^ 
and 9w is the weak mixing angle in the on-shell scheme: 

s 2 w = 1 - M 2 W /M 2 Z (3.4) 

The replacements (3.3) modify the tree level axial vector and vector coupling constants: 



1 

9A = ^ 
9v = 9a 



7ra 
1-4^ 



(3.5) 
(3.6) 



Since however da EC /dc is multiplied, in Eqn.(3.1) by the electromagnetic correction of 
0(a), the effect of the weak corrections in these terms is quite negligible ( < 0.1% ) and 
so da /dc may be considered to be simply Eelectromagnetically Corrected. The label i 
in the last term in Eqn.(3.1) denotes the different partial cross-sections and interference 
terms (see for example Refs. [3,10,30]) resulting from the two (four) Feynman diagrams 
that contribute for / ^ e ( / = e ). Sf BOX are the corrections due to 77 and 7 Z box 
diagrams . 

In (3.1) initial state radiation is exponentiated whereas final state radiation and ini- 
tial/final interference effects are calculated only to 0(a). This simplification is justified 
by the following argument. The large logarithmic terms at 0(a 2 ), 0(a 3 ),... associated 
with initial state radiation do not cancel in the cross-section when the VS and 'hard' con- 
tributions, separated by the arbitrary cut y on the scaled photon energy, are added, and 
give a sizeable (several %) correction. On the other hand, for final state radiation, because 
of the KLN theorem |3l| the leading logarithms cancel exactly in the fully integrated 
cross-section, and cancel almost completely when only loose cuts are applied. For the 
level of accuracy aimed for in BHAGENE3 ( ~ 0.5% fractional error in the cross- section) 
it is then sufficient to consider, in cross-section calculations (both VS and 'hard'), final 
state radiation and initial/final interference effects i effects to O(a). 



The differential cross-section (3.1) is integrated over the angular range cmin < c < 

CMAX '■ 



dc 



5 See Ref. [32] for a discussion of the effect of initial/final interference in the charge asymmetry for 
fermion pair production (/ 7^ e). 
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The probability that a VS type event is generated is proportional to a vs . The lepton 
scattering angle Q\ is generated by first calculating an integrated probability distribution 
in the form of a histogram with content P« in bin % where: 



da vs J ■ 
-de 



/a vs (3.8) 



' c min dc 

Inverting this distribution by linear interpolation yields a LUT for c with bin index j : 

c, = f(Pj) (3.9) 

c is then generated using: 

Cj = f(Rn) (3.10) 

In Eqn(3.10) and subsequently Rn is a random number uniformly distributed in the 
interval < Rn < 1. c is calculated by linear interpolation in the LUT (3.9) using the 
bins k, k + 1, adjacent to P = -Rn (P k < Rn < Pk+i)- 

A similar technique (see also Appendix C) is used throughout the program for the 
generation of one or two dimensional distributions. The advantages are: 

(i) Complete generality. 

(ii) 100% efficiency Multiple trials, as in the weight rejection technique are unneces- 
sary. 

(iii) Fast generation. The simple operations needed to perform linear interpolation in 
a LUT use, in general, much less computing time than the calculation of a weight. 

The main disadvantage is a certain overhead in computing time during the initialisa- 
tion phase, when the LUT's are created. This becomes progressively less important as 
large samples of events are generated. The number of bins used in the LUT depends on 
the steepness of the function to be generated and the desired accuracy. Because of the 
sharp peaking of the angular distribution near c=l in Bhabha scattering, due to t-channel 
photon exchange, a finer binning of the integrated probability distribution (3.8) is used 
in this region. 200 bins are assigned to each of the angular regions: 

cmin < c< 0.8c M ax 

and 

0-8c M ax < c < c M ax 
In the LUT (3.9), 500 uniformly spaced bins are used. 



4 Hard Photon Corrections 



4.1 Generation of Hard Photon Events 



The 'exact' hard photon cross-section is calculated by first generating events according 
to simple factorised differential cross-sections, valid for collinear initial state or final state 
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radiation, and then reweighting each event according to the procedure described, for 
example, in the first of Refs.[4]. 'Exact' is written in quotes because the hard cross 
section formula used is not the result of an exact, fixed order, matrix element calculation, 
but rather an improved soft approximation comparable to that of Ref.[33]. From now on 
this distinction will be understood and the term 'exact' used without quotation marks. 
The exact cross-section formula (see Appendix B) is derived by exponentiating the initial 
state radiation terms of the 0(a) hard cross-section H in such a way that consistency 
with Eqn.(3.1) is obtained in the soft photon limit: 



V 



2E^ 



s < 1 



(4.1) 



The use of such a 'soft' approximation for the generation of 'hard' photons, may be jus- 
tified, in the region of the Z peak, by the small width of the Z relative to its mass : 
Tz/Mz ~ 0.03 which strongly damps out hard initial state radiation. Typically yo in 
Eqn.(3.1) is chosen to be 0.005, corresponding, at the Z peak, to Y^EJ = 225 MeV, 
so that the majority of 'hard' photons are indeed soft at the scale of V z = 2.7 GeV. 
Exponentiation is not applied to the final state radiation terms, since the KLN Theo- 
|3~T| guarantees the smallness of 0(a 2 ), 0(a 3 ),... corrections in this CclSC, clS already 
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discussed in Section 2 above. Similarly, in accordance with Refs. [11-13,32] the initial/final 
interference is also treated only to O(a). 

The approximate differential cross-sections da\ , da^ f° r initial, final state radiation 
respectively, are: 



j asy da (s',t* 

1 167r 3 K + /t_ d{— t*) 



;i + ^)/ e 



, v 

V H 

y 2 



dfl + dQjdy 



da F A 



asy da (s,t) 



dfl+dQ-fdy 



16n 3 k' + k'_ d(-t) 7 
The notation closely follows that of Ref. [3,4]. 4- vectors are defined according to: 



(4.2) 
(4.3) 



and then: 



s = ( p+ +p_) 2 

s' = (q+ + q-f 
k± = p± ■ k 



t= (p+ - q+f 

q± ■ k 



u={p+~ q-f 
u' = (p_ - q+ f 



(4.4) 



dQ + = d(cos9 + )d<p^ 



dQ y = d(cos^ 7 )d0 o 



(4.5) 



6 + , (p + are the polar and azimuthal angles of the Z + relative to the incoming 
direction, while 9 1 , <p 1 are the polar and azimuthal angles of the photon relative to th< 
incoming e + ( / + ) directions for initial (final) state radiation. Also in Eqns. (4.2,4.3): 



—s 



cos e* 



(4.6) 



12 



where 8* is the l + scattering angle relative to the incoming e + direction in the Outgoing 
Dilepton Rest Frame (ODLR). dao is the Born level (uncorrected) differential cross- sec- 
tion for e + e _ — > . The probabilities to generate 'initial state' (I) or 'final state' (F) 
events, according to the approximate models are proportional to the cross-sections a A , 
a A given by integrating Eqns. (4.2,4.3) respectively, over the angles and energies of the 
outgoing leptons and photon: 




vmax dy r c liAx da (s',t*) 
vmin V Jc* MIN d(—t*) 



l + 5V)/ e - y + 



r 



dc* 



(4- 



where 



s(l-y) t* 



y)(i 



cose* 



where 



a 
Ait 



In 



ml 



J Vmax) 



+ 

MIN 



da (s, t) 
d(-t) 



dc^ 



cos 9 i 



(4.9) 



The angular integrals over the lepton scattering angle may be performed analyti- 
cally |0. The results are given below in Appendix A. The expression for the exact 



differential cross-section da EXACT is very lengthy and may be found in Appendix B. 

The approximate cross-section da^, and the exact cross- section da EXACT contain, 
because they are exponentiated, contributions from 2,3,.. hard photons, even though they 
are functions of the angular variables of a single 'photon'. In fact the angular variables 
of all, except one, of the photons are already integrated out. On performing the angular 
integration for this 'photon' i a distribution differential in only the total photon energy 
(the derivative of Eqn.(3.1) with respect to yo, with the replacement yo — > y) is obtained. 
The ansatz employed is to use the same average weight for all pure initial state or all pure 
final state topologies, and, in the case of one initial and one final state photon in the same 
event, to use the harmonic mean of the initial and final state weights. 

The calculation of the average weights W T , [W'j] for events generated according to 
Eqn.(4.2), [(4.3)] proceeds as follows. At the end of the initialisation phase of the program 
N(l, 0),and N(0, 1) events are generated with relative probabilities proportional to a ! A and 
<j a according to Eqns. (4.2) and (4.3). The details of the event generation algorithm used 
in each case are given below in Sections 4.3.1 and 4.4.1. The average weights are then 
given by the expressions: 

i ^(1,0) j^EXACTf^i \ 
W' T = T — ^ ^ (4 10) 



Af(0,l) £ da' A (a{) + da^(ai: 



3 Note that the energy of this 'photon' is in fact the summed energy of all real photons. 
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where a\, {a 3 k ) are the k kinematical variables necessary to completely specify the kine- 
matical configuration of the zth I event (j'th F event). In the subsequent event generation 
phase weights are assigned to multiphoton events according to Table 1. The relative 



probabilities of radiating 1,2,3 photons are given by a Poisson distribution function |35 

pj>27) /P/ (>l7) = 1 _ £ -re ( 4 _ 12 ) 
p£>27)/p(>l7) = 1 _ e -r f (4 _ 13) 

P^/P^ = l-(l + rf ) e - r f (4.14) 

where 

r f = (3 f \n(-) (4.15) 




These probabilities are used to construct the 'a priori' event generation probabilities 
P(rz^,n^) presented in Table 2. The following definitions are used in this Table: 

P vs = a vs /a$ OT (4.16) 

Pi = °a/4ot (4.1T) 
P F = a F Ja* OT (4.18) 



where 

and 

where 



„a vs i / I F 

a TOT = a +cr A + a A 



Pif = (Sj + S F )/{Sj + S F + JS!S F ) (4.19) 



Sj = Pj(l - e~ r *), S F = P F (1 - e~ r f) (4.20) 

In Table 2 P(2, 0), P(0, 3) are actually the probabilities for > 2, (> 3) initial, (final) state 
photons, so so that initial states with 2,3,.. photons and final states with 3,4,.. photons 
are assigned exactly 2 and 3 photons respectively. P(l,l) is derived from P(2,0) and 
P(0, 2) using a factorisation ansatz. The probabilities in Table 2 sum to unity. For the 
case of single photon I or F events the weight W is calculated as: 

W = — [ak) (4.21) 

da A (a k ) + da A (a k ) 

The generation of the kinematical variables a k necessary to completely define the event 
configuration is described below. After checking for consistency with the imposed kine- 
matical cuts (events failing the cuts are assigned weight zero) events with unit weight are 
written out by requiring that: 

W/Wmax > Rn (4.22) 

Wmax is chosen as small as possible (typically Wmax = 2) in order to maximise the 
efficiency of event generation. 



4.2 Cross-Section Calculation 



The calculation of the cross-section printed out after completion of the event generation 
loop is now described. Suppose that N event generation trials are made in the loop with 
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N(nL, n F ) trials in each different final state channel chosen according to the probabilities 
in Table 2. Event weights are defined as follows : 

d EXACT/ i \ 

W i = A if iu^ *V Ki<N{l,0) (4.23) 

J EXACT I 3 \ 

W f = ~r~Tri\ 7pw~Ri Kj<N(0,l) (4.24) 



W^ s = 1, l</<iV(0,0) (4.25) 

where the a]f are defined after Eqns.(4.10,4.11). The cross-section a CUT corresponding 
to the imposed kinematical cuts (trial events that fail the cuts are assigned weight zero) 
is given by: 

a CUT = W(a v > s + a\ + a F A ) (4.26) 

where 



W = [(N(1,0) + N(2,0)W I + (N(0,2) + N(0,2) + N(0,3))W 



F 



+N{l,l)yJW I W F + N(0,0)W vs ]/N (4.27) 



1 JV(1,0) 



W ' = AUO) § W ' (4 ' 28) 



1 JV(0,1) 



W F = — r Y Wi (4.29) 

iV(0,l) fr{ F K J 



| N(0,0) 



Wvs = g < s (4.30) 



The weights W 7 /, W F defined by (3.26), (3.27) are identical to W\, W' f defined by 
(4.10,4.11). Since however, the former are more precise than the latter, due to the typically 
smaller statistical errors of the Monte-Carlo integration in the main generation loop, they 
are preferred for the calculation of a CUT . The error on a CUT , Aa CUT due to the statistical 



error of the Monte-Carlo integration is |36 

Aa CUT 



a CUT 



[SW2 - (SWy/Ni F ] 2 /SW (4.31) 



where 



AT(1,0) N(0,1) 

SW2 = £ (WD 2 + ( W f? 

i=l j=l 

SW = N(l,0)Wi + N(0,l)W F 
N IF = iV(l,0) + iV(0,l) 



15 



4.3 Initial State Radiation Events 



4.3.1 One Initial State Photon 



The final state 4-vectors of events with a single initial state photon are generated by 
the subroutine ZINIGB. These events are first generated according to the approximate 
differential cross-section (4.2) and then re- weighted according to the exact differential 
cross-section using Eqn.(4.23). 

During the initialisation phase a LUT for the scaled photon energy y = E 7 /E, with 
200 uniformly spaced bins from yuiN to Umax is created, as described in Section 2 above. 
This is done by integrating Eqn.(4.2) first over cos 9*, then over y. The angular integral 
is done analytically (see Appendix A), and the y integral by numerical integration @. At 
the same time a two dimensional LUT (200 bins in y, 500 in cos#*) is created. Details of 
the creation and use of such a two dimensional LUT may be found in Appendix C. The 
subroutines used to create these one dimensional (1-D) and two dimensional (2-D) LUT 
are SETYD, SETCTD respectively. 

The kinematical variables used to completely specify the configuration are, in order 
of their generation: 



y, cos# 7 , 0+, cos 9*, 7 

The scaled photon energy y is chosen from the 1-D LUT. The angle (9 7 between the 
incoming e + and the photon is chosen according to the distribution: 

^ - — ^ (4.32) 



where 



c 7 = cos 6* 7 , Pin — \ 1 



Ami 



s 

In this case the integrated probability distribution is calculated analytically and inverted 
to give: 

c 7 = (a - l)/[(3 IN {a + 1)] (4.33) 

where 

a = [(1 - Pin) I {I + /3 IN )\ exp[2Rn\n 1 " >IX 



1-/3 



IN 



The azimuthal angle <p\ between the planes PI (defined by photon and the incoming 
e + ), and P2 (defined by the outgoing Z + and the incoming e + ), in the outgoing dilepton 
rest frame (ODLR), is generated uniformly. Next the 2-D LUT is used, together with the 
previously generated value of y to give cos 9*. Finally the angle 7 , giving the orientation 
of the plane defined by the photon and the incoming e + in the incoming e + , e~ (LAB) 

7 Using the CERN Library program DGAUSS 
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frame, is generated uniformly. To obtain the 4- vectors in the LAB system of the outgoing 
Z + ,/~,7 a Lorentz transformation between the ODLR and the LAB frames is required. 
The corresponding boost direction is parallel to the photon direction. Denoting by as the 
angle between the photon and the incoming e + in the ODLR frame, the transformation 
gives: 



sin as 



cos as 



sin 9 1 
1 — (3* cos 
cos# 7 
1 — (3* cos 



0~ 

PL 

9, 



(4.34) 
(4.35) 



7 



where 

P* = y/(2 - y) (4.36) 

The 3-momentum components of the l + in the ODLR frame: 

g*s*sin</>* (4.37) 
q*s* cos (f>% (4.38) 
q*c* (4.39) 

where 

s* = sin9*, c* = cos6>*, 

are transformed into the LAB frame using the Lorentz transformation defined by Eqns.(4.34- 
4.36). Conservation of 3-momentum in the LAB frame (the 3-momenta of the l + and 
photon being known) then determines the 3-momentum of the l~. Energy conservation 
of the complete event is now checked. In the case of a deviation of more than ±0.5% the 
event is rejected and a new one is generated. A similar check is made for all other hard 
photon event topologies described in the following Sections. If the event is accepted, the 
3-vectors of all particles are rotated about the incoming e + axis, so that the azimuthal 
angle of the photon becomes 7 . 



A 

,(2)* 



Y (3)* 



4.3.2 Two Initial State Photons 



Events with two initial state photons are generated by subroutine ZINI2G according 

14 , [37]] • The differential cross-section is: 



to an improved soft photon approximation 



C n da Q (s') 



P+-P- 



P+-P- 



(P+ ■ fa)(P- ■ h) (p+ • h)(p- ■ k 2 ) 



" E 2E 2>y 



k 2 k\ d?ki d 3 k 2 

~e + wy~k x k~2~ 



4C„ 



dcr (s')d(yi)d(y2)dyidQidy2dQ2 

i - ffi N cos 2 e 7 i)(i - pj N cos 2 e 72 ) 



(4.40) 
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where 2 

s' = M? +l -, E = ^~s/2, y i = 2k i /V^, d{y) = {1 - y + V -)/y 

Here dao(s') is the Born differential cross-section and ki, k 2 ;dfli,dfl 2 are, respectively, 
the energies and solid angle elements of the radiated photons. 9 7 i, 6 l2 are the angles 



between the photons and the incoming e + . The factors d(yi) are the Gribov-Lipatov [p8 | 
kernels. In Eqn.(4.40) recoil effects are neglected, whereas in the event generation al- 
gorithm, the change in direction of the radiator due to the recoil from the first radiated 
photon (in the case that both photons are radiated from the same incoming e*) is allowed 
for. The photon angles relative to the radiator are still however generated according to 
Eqn.(4.40), so that effects due to the virtuality of the first recoiling e are not taken 
into account. Also symmetrisation effects due to different possible time ordering in the 
radiation of two photons from the same incoming line are neglected. 

The kinematical variables used to define the configuration of an event with two 
initial state photons are, in order of generation: 



COS# 7 i, COS# 7 2, M, 7 12, 0+, COS 8 , 



Here y is the scaled total photon energy: 

2y M iN <y = 2(ki + k 2 )l4~s < 2y MAX (4.41) 

u is defined as ln(?/i/?/ 2 ) where y\ and y 2 are the scaled energies of the photons: 

yMiN <yi = 2ki/ y/s < y M Ax (4.42) 

and yi > y 2 . The angles 6> 7 i, l2 are generated according to the same distribution 
(Eqn.(4.32)) as the angle # 7 for the case of single initial state photon events. 7 i2 is 
the angle between the planes defined by the incoming beam direction and each of the 
photons. The angle <p* + is that between the plane formed by the summed momentum 
vector of the two photons and the beam direction, and the plane formed by the Z + and 
the beam direction, all in the ODLR frame. The other variables are defined in the same 
way as for single initial state radiation events. The variables y and u are chosen using 
LUT generated at initialisation. The joint distribution in y and u is derived from that for 
yi and y 2 (Eqn.(4.40)): 



d 2 n 



MAX 



dy\dy 2 

where: 



da {s',c*)dc*d{ yi )d{y 2 ) (4.43) 



s' = s(l -yi-y 2 + ym) ~ s' (4.44) 



In general the outgoing dilepton effective mass v s' depends also on the angles of the 
radiated photons i. The expression (4.44) holds for back-to-back photons, the configura- 
tion corresponding to the maximum value of y in (4.40), and therefore spanning the full 



3 The exact expression, in the present case, is given after Eqn.(4.57) below 
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possible phase space. For collinear photons, where Umax = 1, the yiy 2 term in (4.44) is 
absent. Separate LUT are thus generated for events with (approximately) back-to-back 
or collinear photons. If cos0 7 i, cos# 7 2 have opposite signs (back-to-back configuration) 
then Eqn.4.44 is used to define s' in Eqn.4.44 and Umax — 2. If cos# 7 i , cos# 7 2 have the 
same sign (collinear configuration), the y±y 2 term in Eqn 4.44 is omitted and yuAx — 1- 



By the change of variables: 



Vi = + (4-45) 

y 2 = |(1 -R) (4.46) 

R = ( e u - l)/(e u + 1) (4.47) 

W = ln(l/y) (4.48) 



(4.43) may be re-written as: 



2 H 



d 2 n 



dWdu 

where: 



■MAX , . , . . , . 

da (r' ' • 

MIN 



da (s', c*)dc*d!d 2 (4.49) 



di = 1 - Vi + f , l/i = e u -^/(e u + 1), y 2 = e- w /(e u + 1) 

1-D LUTs for y = e~ w are created at initialisation by subroutines SETYDJ,SETYDL 
for back-to-back, collinear photons respectively using Eqn. (4.49). The integral over c* is 
performed analytically, and those over W and u numerically, using, respectively Simpson's 
rule and Gaussian integration. At the same time 2-D LUTs in the variables W and 
u are created by the subroutines SETRI, SETRL for back-to-back, collinear photons 
respectively. The integration limits on W = ln(l/y) are derived from Eqn. (4.41). Overall 
energy conservation leads to the following limits for the variable u : 

umax = H(y + y ^ AX - y r N )l{y-yf AX + y^ IN )] (4.50) 

umin = -Umax (4.51) 

where: 

y™ AX = mm{y M AX,y - Vmin] 

y^ IN = m&x[y M iN,y ~ yMAx] 

Because the factorised formula (4.49) used to generate the photon energies does not 
take into account the angles of the radiated photons, overall overall energy conservation 
of the event is checked, after generation of the variables a& and before construction of 
4-vectors. In the LAB system, the 3-momentum of the outgoing dilepton is equal to the 
total 3-momentum of the two photons. Energy conservation then requires that: 

yfi = k x + k 2 + [(fci + k 2 f + M, 2 +,_]3 (4.52) 
> h + k 2 + \k\ + k 2 \ (4.53) 
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The inequality (4.53) follows from the condition that the dilepton effective mass M t +i- is 
> 0. (4.53) may also be written as: 

2 > yi + 2/2 + [yl + y\ + 2y 1 y 2 (sin 6 yl sin 6 y2 sin <p yl2 + cos 6> 7l 6» 72 ] * (4.54) 

If the R.H.S. of (4.54) is found to be > 1.99 the current event is rejected and a new one 
generated. 

In constructing the 4-vectors of the photons, the recoil of the radiating e ± is allowed 
for. The recoil angle, a R of the e ± , after photon 1 has been radiated, is: 

/ yi sin d-yi \ 

a R = arctan : — — r (4.55) 

\ 1 - y 1 1 cos yl | J 

In the case of radiation of the two photons from the same incoming line, given by the 
condition: 



cos 7 i cos # 7 2 > 

the angle 6 l2 is chosen relative to the direction of the recoiling e ± , rather than the 
beam direction. 

The remainder of the event generation follows closely the procedure for the case of one 
initial state photon, described above. The lepton scattering angle is chosen using the same 
2-D LUT as in the one photon case. The boost between th LAB and the ODLR frames 
is now along the direction of the vector sum (K) of the 3-momenta of the two photons. 
The angle 6 1 in Eqns. (4.34,4.35) is replaced by 9k, the angle between the incoming e + 
direction and K. The Lorentz transformation parameter f3* of Eqn.(4.36) is replaced by 
the expression: 

^ ^ vH^) (4 ' 56) 

and the lepton momentum in the ODLR frame used in Eqns. (4.37-4.39) is now given by: 



q* = fj-ml (4.57) 

where 

s' = - A {2-yf-{Kf 

Finally, since events were generated using the approximate value of s', s' given by Eqns. (4.43) 
and (4.44) (with the appropriate modification for collinear photons) the events are re- 
weighted using a WRP, to take into account the exact value of s' as given above after 
Eqn.(4.57). The weight used is: 

da (s',c*)/da (s',c*) 
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4.4 Final State Radiation Events 



4.4.1 One Final State Photon 

Events with a single final state photon are generated by the subroutine ZFINGB, 
according to the differential cross-section Eqn.(4.3). They are subsequently re- weighted 
according to the exact differential cross-section using Eqn.(4.24) 

The kinematical variables a:*; used to define the event configuration are, in order of 
their generation: 

y, cos 9*, cos9±, 0^, <f>± 

Following (4.8) the scaled photon energy, y, is given by: 

V = yMiNexp[Rnln I Vmax ]] (4.58) 
V Vmin ) 

The angle, 6*, between the / + and the photon in the ODLR frame, is chosen according to 
the distribution: 

*L 1 U 59) 

where 



c* = cos6>* (3_ 



FI 



2 
S 



The scaled energies of the / + , / in the LAB frame, y+, y_, are obtained by a Lorentz 
transformation: 

y ± = ^iq^±[3*q*c;] (4.60) 
V s 

where 



2 



ir = v/(2-v), 7* = (i - fVV 1 - y 

The LAB scattering angles 9± of the are found using the same 2-D LUT in y and 
cos#* as used for initial state radiation. Here y is set to zero so that s' — > s and 9* — > # + . 
Following Ref[14], the angle #+ is assigned to the Z + , or 7r — 9 + to the Z~, according to 
whether, or not : 

P s > Rn (4.61) 

where 

Ps = yl/(yl + yl) 
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This has the effect that, if the photon is radiated from the then the scattering angle 
is most likely assigned to the other lepton Z T , whose direction is unaffected by recoil 
effects. In the case that the 9 + is assigned to the / + , the 3-axis is chosen along the l + 
direction and the 3-momentum components of the / + ,7 are given, in the LAB system, by 
the expressions: 

= 

= (4.62) 



(3) 

Q+ 



k^ 



/ xf 

— ys 1 cos0 7 
— ys 1 sin 7 



(4.63) 



where 



< = 1 - 2[1 - (1 - y)/y + ]/y, < = - (d,) 2 , ft = 



4mf 



Here 7 is the azimuthal angle of the photon about the / + direction. If 7r — 9 + is 
assigned to the l~, the 3-axis is chosen along the l~ direction, and the 3-momentum 
components of the l~ , 7 are given in the LAB system as: 






n 



(4.64) 



k^ 



// xi 

— ys 1 cos0 7 

// • xi 

— ys 1 sin 7 



(4.65) 



where 



< = l-2[l-(l-y)/y_]/y, 



sy 2 _ 



The momentum components in the standard LAB system with 3-axis along the in- 
coming e + direction, are found by rotating the above 3-vectors by 9 + (or n — 6 + ) about 
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the 1 -clXIS, clS appropriate. The azimuthal angle </> + ,(0_) of the Z + , (/ ) about the 3-axis 
is then generated uniformly. Finally the 3-momentum of the remaining lepton is given by 
momentum conservation in the LAB frame. 



4.4.2 Two Final State Photons 

Events with two final state photons are generated by the subroutine ZFIN2G according 
to the factorised differential cross- section §, |l|, ^7|: 



j f n A is Q+-Q- 
da 2ry = C n da (s) 



n UV ' (<?+ ■ fci)(?_ • h) (q + ■ fc 2 )(g_ • k 



2 



k\ k\ k 2 k\ d 3 kid 3 k 2 

X{1 ~E + 2^ )(1 " E + 2&> h k 2 

da (s)d(y 1 )d(y 2 )dy 1 dn 1 dy 2 dn 2 

° n (i -p FI cos 2 e;o (i-p 2 FI cos 2 e; 2 ) 1 m} 

The variables in (4.66) are defined after Eqns.(4.40) and (4.59), with the exception of 9' x , 
6 ' 2 , which are the angles between the photons 1,2 and the Z + in the LAB system. To take 
into account recoil effects however, when cos^ 2 < tli e photon 2 is assigned the angle 
7i — 6' l2 relative to the l~ direction rather than 9' 2 relative to the / + direction. 

The kinematical variables in the order of generation are: 



y, cos9' 7l , cos9' l2 , u, cos9±, 0^ 12 , <p' lT , <f>± 

Here 9± are the scattering angles of the relative to the incoming e + in the LAB 
frame. The azimuthal angle ft 12 is chosen to be that between the planes defined by : 

(/c 2 ,g+) for cos^ 2 > an d between the planes: (ki,ql); (/c 2 ,gl) for cos^ 2 < 0- 
ft lT is the azimuthal angle about the l + direction of the total momentum vector K 2 of the 
two photons. 4>± are the azimuthal angles of the I relative to the incoming e + . All three 
angles are generated uniformly in the interval to 2n radians. The polar angles 6^, & l2 
are generated according to Eqn.(4.32) with the replacements 6> 7 — > 9' , Pin — > Pfi- The 
variables y and u (y — y\ + y 2 , u — ln(yi/y 2 )) are chosen from, respectively, 1-D and 2-D 
LUT's generated at initialisation. The distribution (4.49) is here replaced by: 

did 2 (4.67) 



dWdu 



where W = ln(l/y). The subroutines that generate the 1-D LUT for y and the 2-D LUT 
for W are SETYDK, SETRK respectively. The angle 9 + is generated as described above 
for single final state photon events. 

The construction of the 3-vectors of the final state particles is now described. Unlike 
for the case of single final state photon radiation it is here more convenient to work entirely 
in the LAB frame. Choosing the 3-axis along the Z + direction, with the 2-axis in the plane 



23 



containing g+ and k\ the momentum components of photon 1 are: 



k? = 





k? = 


^-l/i^i (4-68) 


k\ 3) = 


2 yi^i 


where 






c 7l = cos 6» 7l 


In the case that c' 72 > (photon 2 radiated from the / + ) the momentum components 
of photon 2 are given by: 


fc 2 - 2 


y 2 s 7 2 cos 7l2 


1.(2) 

^2 - 2 


?/ 2 s 72 sin 7l2 (4.69) 


1.(3) _ 

2 ~ IT 





For c 72 < (photon 2 radiated from the l~) the direction of the 3-axis in Eqn.(4.69) is 
chosen opposite to the direction of the l~. In this case the momentum components are still 
constructed according to Eqn.(4.69), but they are subsequently rotated about the 1-axis 
through the acollinearity angle a a between the l + , l~ directions. This angle is calculated 
from the equations: 

sin a A = yis'^/y- 

y. = 2-y 1 -(l-y 1 )/[l-0.5(l-c / 7l )y 1 ] (4.70) 



cos a a = y 1 — sin 2 a a 

By this procedure the strong peaking of the radiated photons along the directions is 
properly accounted for, even in the case of recoils generated by the first radiated photon i. 
Before proceeding to construct the 3-vectors of the outgoing leptons energy conservation 
is checked using Eqn.(3.51). If the R.H.S. is > 0.99^ the event is rejected and a new 
one is generated. 

The scaled energies y± of the are next calculated: 

V + = 2[l-,(l-f)-^£]/(2-, + ^) (4.71) 
4 s V s 

? y_ = 2-y-y + (4.72) 



9 In Eqn.(4.69) the angle 0' l2 is that between the photon and the radiating (virtual) l~, rather than 
that between the photon and the outgoing l~ as in Eqn.(4.66). This difference is of little importance for 
soft or nearly collinear photons. For hard non-collinear photons however the rate will be over-estimated 
as compared to Eqn.(4.66). In the case that both photons are radiated from the l + the ansatz used tends, 
on the contrary to underestimate the rate of such photons as compared to Eqn.(4.66) 
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As in Eqns. (4.68,4.69) the 3-axis is along the / + direction. The angle 9+ is assigned to l + or 
7r — 9 + to the l~ , as described above for single photon final state radiation. The azimuthal 
angle of the total photon momentum K 2 is then generated uniformly by assigning the 
angle <f)' lT . The 3-axis is now rotated into the direction of the incoming e + and the angle 
(f) + (or 0_) is assigned. Finally the 3-momentum of the remaining lepton is found using 
momentum conservation in the LAB system. 

4.4.3 Three Final State Photons 

Events with three final state photons are generated, according to the obvious gener- 
alisation of Eqn.(4.66), by the subroutine ZFIN3G. The event configuration is defined by 
the following kinematical variables: 

V, cos 6^, cos cos 9' y3 , it', u, cos6> ± , 0' 12 , 0^, (f)' jT , <p± 

6^3 is similarly defined to 9' 1: 9' j2 but refers to the third (least energetic) photon. The 
variable u' is given by: 

The definition of is similar to that given above for the case of 2 final state photons 

— * 

except that the total photon momentum vector K 3 corresponds to three photons. The 
angle is that between the plane defined by the total momentum of the first two (most 

energetic) photons K 2 and the / + (or / _ ), and that defined by the direction of photon 3 
and / + (or l + is taken for cos 6^ 3 > 0(< 0). All other variables have been defined 
previously. 

Following the soft photon limit of the generalisation of Eqn.(4.66) the photon energies 
are first generated according to the distribution: 



d 3 n 



(4.74) 



dyidy 2 dy 3 ymyz 
By a change of variables: 

W = ln[l/(yi + y 2 + y 3 )], u = \n(y 1 /y 2 ), u = ln[(j/i + y 2 )/y 3 ] 
(4.74) simplifies to 

d 3 n tA 

— - — ~ constant (4.75) 

dWdudu' 

A 1-D LUT for the variable W, and a 2-D LUT for vl (for a given W) are created 
at initialisation by the subroutines SETYD3, SETRK3 respectively. The limits on the 
variable u in the nested integration: 



dW i / du' 



du 
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are 



(y> _ y MAX\ 

u MIN = In ( y M l AX j (4.76) 
Umax = In ^ y ^ N j (4.77) 



where 



yf /7V = max[?/' - i/jvmi, 2/m/at] 

y MAX _ min ^/ _ y M/JV) y MAX ] 

and 

y' = 2/i+2/2 = 2/e u 7(l + e u ') 

The limits on the u' integration are given by Eqns. (4.50,4.51) with the replacements 
u —>■ u' and 

y MiN _^ max [2y M/7V) y - y MAX ] 
y MAX _^ m i n [y MAXj y - 2y M w] 

The integral over u is done analytically, that over u' by Gaussian integration, and that 
over W by Simpson's rule. During the event generation phase the value of W is chosen 
from the 1-D LUT, The value of u' is then chosen using the 2-D LUT. Finally, u is given 
by: 

u = umin + Rn(u M iN - Umax) (4.78) 

where umin, Umax are given by Eqns. (4. 76, 4. 77). yi, y2, y% are then derived from W,u, 
u' via the equations: 

w 

y = e 

y' = ye u '/(l + e u ') 

Vl = y'e u /(l + e u ) (4.79) 

V2 = y'-y 
ys = y-y' 

The photon spectra are then modified by weight rejection to take into account the hard 



photon corrections given by the Gribov-Lipatov [5^] kernels. 
Defining the weight function: 

W = d(y 1 )d(y 2 )d(y 3 )y 1 y 2 y 3 (4.80) 
where d(y) is defined after Eqn.(4.40), the event is rejected if: 

W < Rn 

The construction of the momentum vectors of the final state particles follows closely that 
described in the previous section for the case of two photons. The 3-vectors of the two 
most energetic photons are given by Eqns(4.68,4.69) above, and the angle 6' l2 is chosen 
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relative to the Z + direction (opposite of the l~ direction ) according to whether cos 9' l2 is 
> or (< 0). Energy conservation is then checked using Eqn.(4.54). To take into account 
recoil effects in the radiation of the third photon the acollinearity angle between the l + 
and l~ after the first two photons have been radiated is calculated, and 9' l3 is chosen 
relative to the Z + direction (opposite of the l~ direction, for cos9' l3 > (< 0). That is, 
the procedure used above for the second photon is iterated. The azimuthal angle 73 is 
chosen to be that between the planes defined by (K 2 ,q+) ; (£3,9+) for cos# 73 > 0, and 

that between (K 2 ,ql) ; (k 3 ,ql) for cos 9' y3 < 0. Since the relative directions of all three 
photons are now fixed energy conservation is again checked using the generalisation of 
Eqn.(4.53): 

v^> h + k 2 + h + \h + k 2 + k 3 \ (4.81) 
If the RHS of Eqn(4.81) is > 0.99^, the event is rejected. Eqns. (4.71,4.72) with the 

— * — * 

replacement K 2 — > K 3 is now used to calculate y + , y_. The scattering angle 9 + is assigned 
to the / + , or n — 9 + to the l~, according to Eqns. (4.61). The angles (p' lT and 0+ or 0_ 
are assigned as described above for the two photon case. Finally the 3-momentum of the 
remaining lepton is calculated from overall momentum conservation. 



4.4.4 Events with One Initial State and One Final State Photon 



Events with one initial state and one final state photon are generated by the subroutine 
ZINF2G according to the differential cross-section (c.f Eqns. (4. 40,4. 66)): 



da I F 2l = C n da (s") 



P+-P- Q+ ■ Q- 

(P+ • fci)(P- • fci) (q+ ■ k 2 )(q- ■ k 2 ) 
h kf k 2 k\ <Pki<Pk 2 

, r da (s")d(y 1 )d(y2)dyidn 1 dy 2 dn 2 

n (i-^cos2e 7l )(i-^ /C os20; 2 ) l ' ] 



where 



s" = s(l- yi ) 

and photons 1,(2) are radiated in the initial, (final) state. The kinematical variables used, 
in this case, to define the event configuration are: 



y, cos9 jl: cos9' l2 , u, cosQ'[, 0", 0" 2 , 7 i 

Here 9" 2 is the angle between photon 2 and the / + in the rest frame of l + l~^2 (ODLGR 
frame). In this frame the 1-axis is chosen perpendicular to the plane defined by the 
incoming e + and the direction of photon 1. 9'\_ , <f>+ are the polar and azimuthal angles of 
the / + relative to the incoming e + direction (allowing, if necessary, for the recoil generated 
by photon 1) in the ODLGR frame. The other variables have been previously defined. 

The variables y, u are generated using 1-D, 2-D LUT created at initialisation by the 
subroutines SETYDI, SETR respectively. The procedure is the same as that described in 
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Section 4.3.2 above. The 2-D LUT is created according to Eqn.(4.49) with the replace- 
ment: s' — > s". The configuration of / + /~72 is generated as described in Section 4.4.1, for 
single final state radiation, but in the ODLGR frame rather than the LAB frame. Recoil 
effects from the initial state photon are accounted for by rotation by the angle 0.3 given 
by Eqns. (4.34-4.35) with the replacement y — > y\ in the formula (4.36) for (3*. Since the 
scaled photon energies y%, y 2 are defined in the LAB frame, the energy of photon 2 must 
first be calculated, by Lorentz transformation, in the ODLGR frame, from its angles in 
this frame (specified by 9'' 2 , 4>'l/2-> 0+ ) an d its LAB energy ^2/2 ■ For simplicity, in 
this case, Q'\ is always assigned to the l + , rather than Q" + to the l + or 7r — d'\ to the l~ 
according to Eqn(4.61) 0. Finally the 4-vectors of l + , l~ and photon 2 are tranformed 
back into the lab. frame and the event is rotated about the incoming e + direction so that 
photon 1 has azimuthal angle <p 7 \. 



5 Program Structure and Performance 



5.1 General Organisation. How to use the Program 

The program has a very short main program containing definitions of the most im- 
portant input parameters, which are stored in the labelled common block ICOM. These 
variables are described in Table 3. The execution of the program has three distinct phases: 

• Initialisation 

• Generation of a single unit weight event 

• Termination 



Each of these phases is entered via a call to subroutine BHAGENE3 in the main program: 

CALL BHAGENE3(MODE,CTP1,CTP2,CTM1,CTM2,CTAC,EPO,EMO) 

MODE is set to —1,0, 1 for the initialisation, generation and termination phases respec- 
tively. The other parameters of BHAGENE3 define the kinematical cuts to be applied to 
the generated events: 

CTP1 = minimum value of cos 61+ 
CTP2 = maximum value of cos 61+ 
CTM1 = minimum value of cos Q\- 

10 The procedures described in Ref.[14] for assigning the scattering angle in the case of single initial or 
final state photons are not directly applicable in this case 
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CTM2 = maximum value of cos 61- 
CTAC = maximum value of cos (p co i 
EPO = minimum energy of l + (GeV) 
EMO = minimum energy of l~ (GeV) 



All these cuts are applied in the laboratory (incoming e + , e~ centre of mass) system. The 
angle 4> co i is the collinearity angle between the / + and the l~ (CATC = -1 for a back-to- 
back configuration). In the calls of BHAGENE3 with MODE = 0,1 only this parameter 
need be specified. A typical main program to generate 5000 unit weight events might be: 



PROGRAM BHAMAIN 
IMPLICIT REAL*8 (A-H,0-Z)) 

COMMON/ICOM/OIMZ,OIMT,OIMH,OMAS,IOCHA,IOEXP,OW,OCTCl,OCTC2,IOXI 



OIMZ=91.18D0 

OIMT=150.0D0 

OIMH=100.0D0 

OMAS=0.12D0 

IOCHA=l 

IOEXP=l 

OW=91.00D0 

OCTC1=-0.8DO 

OCTC2=0.8D0 

IOEXI=713883717 



Z Mass (GeV) 

Top quark mass (GeV) 

Higgs boson mass (GeV) 

Alphas 

for Muon pairs, 1 for electron pairs 
for O(Alpha), 1 for exponentiation 

Collision energy (GeV) 

Lower cos(theta) limit in ODLR frame 

Upper cos(theta) limit in ODLR frame 
! Intial random number 



KINEMATIC AL CUTS FOR GENERATED EVENTS 



CTP1=-0.7D0 

CTP2=0.7D0 

CTM1=-0.7D0 

CTM2=0.7D0 

CTAC=-0.9D0 

EP0=2.0D0 

EM0=2.0D0 



Lower cos(theta) for 1+ in the lab frame 
Upper cos(theta) for 1+ in the lab frame 
Lower cos(theta) for I- in the lab frame 
Upper cos(theta) for I- in the lab frame 
Upper cosine of collinearity angle 
Minimum energy of 1+ (GeV) 
Minimum energy of I- (GeV) 



INITIALISATION PHASE 



CALL BHAGENE3(-1,CTP1,CTP2,CTM1,CTM2,CTAC,EP0,EM0) 



29 



! GENERATION LOOP 



DO J=l,5000 

CALL BHAGENE3(0) 

ENDDO 

! TERMINATION PHASE 

CALL BHAGENE3(1) 

STOP 

END 

Other initialisation parameters of interest to users are defined in BHAGENE3 itself. A 
list of the most important of these can be found in Table 4. 



5.2 Initialisation Phase 

A flow chart of the initialisation phase of the program is shown in Fig I. The func- 
tionality of the different subprograms shown there has been described above. The main 
physical quantities calculated are also indicated. In order to estimate the average weights 
of events with > two hard photons the initialisation phase actually includes the generation 
of 4 x 10 4 events, which is relatively time consuming. Users should be aware of this. 

5.3 Generation Phase 

As shown in the flow chart in Fig 2. each of the different event topologies is generated 
by a different subprogram. The calculations performed by the five different hard photon 
subgenerators : ZINIGB, ZFINGB, ZFIN2G, ZFIN3G, and ZINF2G are described above 
in Section 4. The 4-vectors of the generated events are written in the labelled common 
block C4VEC : 

COMMON/C4VEC/PPV(4),PMV(4), QPV(4),QMV(4),GAM1V(4),GAM2V(4),GAM3V(4), 
WEIGHT,NPHOT,ISG 

The 4-vectors (defined as: p x ,p y ,p z , E ) are in the order: (e + , e~, / + , l~, 71, 72, 73 ). 
WEIGHT is the event weight, NPHOT the number of filled photon 4-vectors and ISG a 
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flag indicating the subgenerator that produced the event (see Table 2.). The photons are 
ordered in energy, the first photon being the most energetic. 

5.4 Termination Phase 

A flow chart of the termination phase is shown in Fig 3. The exact cross-section (a CUT ) 
and its error (Aa CUT ) are printed out, together with the input parameters (including those 
calculated in the initialisation phase). Other cross-sections used to calculate the event 
generation probabilities P(n*n^) are also printed out. A sample output is shown in Fig 
4. 



5.5 Program Performance 

As mentioned above, after a relatively lengthy initialisation procedure, during which all 
look-up tables are created and average weights are calculated for multiphoton events, the 
event generation procedure is itself fast. Typical times for the initialisation phase are 140, 
48 IBM 3090 CPU seconds for e + e~, fi + n~ pair generation, respectively. Average times 
to generate a single unit weight e + e~ , event are 1.41 x 10~ 3 , (0.64 x 10~ 3 ) IBM 

3090 CPU seconds. The combined weight distribution for initial and final state radiation 
e + e~ events for parameters and cuts as in the sample main program given above, (whose 
output is shown in Fig. 4) is presented in Fig. 5. The weight distribution can be seen to 
be well centered around 1.0, resulting in an efficient generation of unit weight events by 
the weight throwing procedure 0. The fractions of events with weights > 2.0, (3.0) is 
~ 6 x 10" 3 , (2 x 10^ 4 ) respectively. The maximum weight may then be chosen to be 2.0, 
allowing efficient generation of event samples with cross-sections known at, or below the 
% level 0. 
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Appendix A 

Following Ref.[39], the Born differential cross-section for Bhabha scattering, including 
both s and t channel Z exchange may be written as: 



^ = ^[5o + 5 2 (-) 2 + 5 3 (l + -) 2 ] (Al) 
dt s z s s 



where 



t = _£(l-cos0+) 

Bo = [- t +c^ X t(t)} 2 
B 2 = \l + c_ Xs (s)\ 2 

B 3 = i{|l + j + a + [ X t{t) + Xs(s)]\ 2 + {a+ -> a_)} 



c- = g v - g A , a ± = (fiv ± #a) 

5 S 
Xt(t) = TJ2 , Xs{s) 



t-M 2 z ' ASW s - Mf + is{V z /M z ) 

9 + is the scattering angle between the incoming and outgoing e . Mz, Tz are the 
mass and width of the Z. The Z width is neglected in the t-channel exchange amplitude. 
The vector and axial- vector coupling constants g A , gv are given, in the Standard Model 
by Eqns. 3.3-3.6. The cross-section for /x-pair production is recovered on setting s/t to 
zero in B , B 3 . 



Analytical integration of Al over t yields the result |31] ( see also the first of Ref. [4]): 
i-tMAx _ 2-na 2 

' IN 



/ da = — H-STW) - S(t MAX )] (A2) 



The function S(t) is the sum of the following 5 terms derived from Eqn. Al : 



, 1 2c_ t-Mf c 



jB dt = -A-j + ^n 



2 



z t t-M 2 z 



1 r 2r 2 2(a+ + a_), t - Ml (a 2 ,+a 2 

-is 1 ^— ^ n — — — -± 

2 1 1 t Mf t t-M 2 

2< 



+2s[{b + + 6_) hit + (a + 6 + + a_6_) ln(t - Mf)] + + b'_)t} 
f^-dt = s{2 hit + 2(a + + o_) ln(t - Mf) + (a^ + a 2 _)[ln(t - M 2 ) - M f ]} 



t 2 

+2{{b + + &_)* + {a + b + + a_6_)[t + Mf ln(t - M 2 Z }} + + b'_) — 



r t 2 
jB 2 -dt 



t 2 B^_ 
3s 2 
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t 2 

B 3 -dt 

s 2 



-{[2t + 2(a+ + a-)(t + M| ln(t - M 2 Z )) + (a 2 + + a 2 _)(t + 2M| ln(t - M%) 

M 4 2 t 2 t 2 

-— )] + -[(&+ + &_)- + (a+6+ + a_6_)(- + tAf| + M\ ln(f - Mf ))] 

z 



where 



6± = 1 + a±ReXs 

b' ± = 1 + a ± 2 J Rex s + a 2 ± \x s 



Appendix B 

The exact hard photon cross-section, with exponentiated initial state radiation, is 
derived from formulae given in the second of Refs. [4] : 



da EXACT a 3 y 



dQ+dQjdy 1Qtt 2 s 



X(a k ) 



(Bl) 



The function X(ak) of the kinematical variables a k specifying the event configuration, 
(defined in Sections 4.3.1 to 4.4.4), has contributions from s-channel exchanges (ss), 
t-channel exchanges (tt) and s — t interference (st). In each of these cases separate 
contributions from initial state radiation (INI), final state radiation (FIN), and initial/final 
interference (INT) may be distinguished. Thus : 



X = Xi i = ss, tt, st; j = INI, FIN, INT (B2) 



where 



X 



INI 



X 



FIN 



X 



INT 



s'k + k_ 
-m 2 J(s 

1 



fli(c_, s')[f{s)t 2 + t' 2 } + B 5 (s')[f(s)u 2 + u 



B ss (s',t,u) B ss (s',t',u') 



/2l 



K 2 _ 



K. 



+ 



B 1 (c_ , s) [t 2 + t' 2 } + B 5 (s) K + u 



/2l 



— m 



B ss (s,t,u') B ss (s,t',u) 

/o ~r /o 



SS 



U U 
+ ■ 



t 



t' 



(,- (,- ^ (,- I/"^ (,- I, - ^ I,- ly'f 

ri/_(_ rv rb — rv_|_ nj_|_nj_|_ rv — rv 
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x 'B 4 (s, s')(t 2 + t' 2 ) + B 6 (s, s')(u 2 + u' 2 
P+ ■ (q+ x Q-)(s ~ s') D /w 2 /2 



2^/is' K+K_ «'+«'_ y 



„, / , r /(sW 2 + U /2 1 „. /Xr /(S)S 2 + S' 2 . 

B e (t,t')[^— ) ]+S 4 (*,0[- 



-m 2 J{s) 



tt' 

B tt {s',t,u) B u (s',f,u') 



Av_i_ AC 



-m 2 



B 6 (t,t')[^^]+B 4 (t,t')[^-} 



B tt (s,t,u') B tt (s,tf,u) 
"i ' "i 



s 2 + s' 2 

tt' 

u 2 + u 

+— 

■u 



._*- r s 1 ( c _,0--^Vs (C _,o 

AC_|_AC_|_ AC — AC 



+ 



+ 



U' 



K-K. 



AC_^_/v_ AC_ AC_j_ 



5a(M0[ !i 1 ^]+54(M')[^ 1 ^] 



X. 



/AT/ 



X 



FIN 

st 



f{s)u 2 + u' [ 



S'K+K- 

-m 2 J(s) 



u 2 + u 12 

SKi\ AC 



-m 2 f(s) 



u' + t' 
t' 



B 6 (s',t') 



u + t 



B 6 (s',t) 



B st (s',t,u) B st (s',t',u') 



AC_ 



B 6 (s,t)-(^±^j B 6 (s,t') 



B st (s,t,u') B st (s,t',u) 

fn ~T fn 



u' 



AC AC I 



+ 



v! + s' 



AC_|_ AC_j_ 



'u 2 + u' T 



+ 
+ 
+ 



u u + s' 

— r + r 

AC_|_AC AC — AC 

■u' u' + s 



AC_|_ AC_ 



It 



ACj_ AC_ 



+ 



■u + s 

AC j_ AC_l 



s't' 
'u 2 + u' 2 " 

v , 

'u 2 + u' 2 ^ 

vw 2N 

St' , 



B 6 (s',t) 



B 6 (s,t) 



B 6 (s,t') 



i ^[pl-(q + xql)](u 2 -u' 2 ) [ ft-t\ BiM | o 



2k + 



tt' 



St(K_K' + K'_) 



+2 



B 7 (s,t') 



-2 



- 2- 



St' (k + k' + k'^_) s't(K + K-K,'_) s't'(K + K-K' + ) 
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£i(c,s) 
B 4 (s, s') 



B 5 (s) 



l + l-(s-M 2 z ) + c 2 ]\xs(s)\ 2 
s 

i+ c -[(i-^)ix s ( S )r+(i-^)ix s ( S ')i 2 ] 

+c 2 .[(i - %{i - ^f) + ^]|x s ( S )l 2 lx s (*')l 2 

M 2 

1 + [2(<£ + ^)(1 - -*) + + ^ + 6<&£)]|x,(a) 



56(5,/) = l + (^ + ^) 



B ss (s,t,u) 
B tt (s,t,u) 

B st (s,t,u) 



{l-^)\Xs{s)? + {l + ^f)\Xs{s')? 
s s' 



= +(9v + 9\ + Q9v9 2 a) 



Ml.. Ml. M 2 T 2 
1 -)(1 + Z , Z 



S 7 (s,s') = -AM z r z \g Agv [ 



s s 

12 |„ ^ /M2 



|X S («)| 2 | Xs (s')P 



s' 



+2^(^ + ^)(i-^)|x s ( S )| 2 |x,( S ')| 2 } 



B 1 {a 2 + ,s) + B l {a 2 _,s) 



[y + 2B 1 (c_,s)( t -) 2 



B 1 (a 2 + ,t)+B 1 (a 2 _,t)} (-) 2 + 2B 1 (c_, t)( 6 -) 2 



u 



45 6 (s,t)- 



\Xs(s)\ 2 \Xs(s- 



l\\2 



The constants c_, a± are defined in Appendix A and g v , g A in Eqns 2.3-2.6. Note that 
the functions Xs, Xt are assigned according to the arguments of B 1 , B 4 ... . For example 
B±{s,t') is given by the replacement Xs(s') —> Xt(t') in the expresssion above for B^(s, s'). 

Exponentiation of the initial state radiation is included via the function: 

/(s)=2C^exp&ln(y)-l (£3) 

where /3 e , C v are defined after Eqn 2.1. This ensures that the y — > limit of Eqn. Bl is 
identical to the derivative of Eqn 2.1 with respect to yo, in which the replacement yo — * V 
is made. That is 'soft' and 'hard' photons are treated in a consistent way. 

Appendix C 

Two dimensional look up tables are generated according to a generalisation of the 
procedure described in the text (Eqns 3.7-3.10) for a one dimensional look up table. 
Consider for example the case of the Born differential cross-section da /dc for a range of 
different CM energies : s m i n < s < s max . Let k be the bin index for s. The equations 2.7, 
2.8 generalise to: 



On = 



dc 



°-dc 



(CI) 
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where 



dal = da (s k ) 



<^0 Jc min "l- 



For each value of k the distribution Pf is inverted by linear interpolation to yield a look 
up table with bin index j : 

4 = /*(#) (C3) 

To generate the angular distribution for a given value of s, the adjacent bins in s of 
index k, k + 1 are first located : 

s k < s < s k +i 

The closest bins in the look up tables of index k, k + 1 to the random number Rn are the 
found: 

P t k <Rn< P? +1 
P k+1 <Rn< P$Xl 

With: 

S t = (Rn - P?)/(P? +1 - P k ) 

5, = (Rn-P^)/(P^-P^ +1 ) 
Two values of c are now calculated according to : 

c k = f k (p t k ) + im+i) - f k (p*)\5i 

c k+1 = f k+ \p k ,) + [f k+1 (P^D - f k+1 (PZ +1 )}Si> 

with 

A fe = (s - s k )/(s k+1 - s k ) 
the generated value of c is, finally, given by the equation : 

c = c k + (c k+1 - c h )A k (C4) 
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Assigned Weight 


2 





W'j 





2 


w' F 





3 




1 


1 





Table 1: Weights assigned to multiphoton events 







P( < , < ) 


ISG 








Pv,s 





1 





Pie** 


1 





1 


P F e-' r f 


2 


2 





- e~ r <)p IF 


3 





2 


P F r f e- r f)p IF 


4 





3 


P F [l-(l + r f )e- r f]p IF 


6 


1 


1 


[PjP F (l-e- r «)(l-e- r O]W 


5 



Table 2: 'A priori' probabilities for different hard photon multiplicities 



OIMZ 


Z mass (GeV) 


OIMT 


Top quark mass (GeV) 


OIMH 


Higgs boson mass (GeV) 


OMAS 


a s (M z ) 


IOCH 


= 0(/i+/!~),= i(e + e-) 


IOEXP 


= f exponentiated , = 0(a) 


OW 


collision energy (GeV) 


OCTCf 


lower cos 9i+ in the ODLR frame 


OCTC2 


lower cos Q\- in the ODLR frame 


IOXI 


initial random number 



Table 3: Variables of the labelled common ICOM. 0CTC1,0CTC2 are used in setting 
up the LUT of the lepton scattering angles. To allow for the effects of the Lorentz boost 
the angular range should be chosen somewhat wider than that defined by the cuts in the 
LAB system. 
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NPAR(1 
NPAR(2 
NPAR(3 
NPAR(4 
NPAR(6 
XPAR(1 
XPAR(2 
XPAR(3 
XPAR(4 
XPAR(9 
XPAR(IO) 

YMA 

YMI 
WTMAX 



1, weak loop corrections ON, OFF 
2,3 parameterisations of had. vac. pol. 
0,1,2 two-loop aa s ml correction 
1,0 weak box diagrams ON, OFF 
1, two- loop terms oc raj ON, OFF 
initial lepton charge (-1.D0) 
final lepton charge (-1.D0) 
final lepton colour (1.D0) 
final lepton mass (GeV) 
QCD correction to (non-6 quarks) 
QCD correction to Tf 
maximum value of J2E-y/Eb eam (0.99D0) 
minimum value of E 7 /E beam (0.005D0) 
maximum value of the event weight (2.0D) 



Table 4: Control parameters defined in SUBROUTINE BHAGENE3. Default values are 
underlined or given in parentheses. 
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FIGURE CAPTIONS 



Fig 1 Flow Chart of the Initialisation Phase. 

Fig 2 Flow Chart of the Event Generation Phase. 

Fig 3 Flow Chart of the Termination Phase. 

Fig 4 A typical line-printer output from BHAGENE3. 

Fig 5 Distribution of weights for events with hard initial or final state photons. 
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